An efficient method to calculate excitation energy transfer in light harvesting 
systems. Application to the FMO complex. 
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A master equation, derived from the non-Markovian quantum state diffusion (NMQSD), is used to 
calculate excitation energy transfer in the photosynthetic Fenna-Matthews-Olson (FMO) pigment- 
protein complex at various temperatures. This approach allows us to treat spectral densities that 
contain explicitly the coupling to internal vibrational modes of the chromophores. Moreover, the 
method is very efficient, with the result that the transfer dynamics can be calculated within about 
one minute on a standard PC, making systematic investigations w.r.t. parameter variations tractable. 
After demonstrating that our approach is able to reproduce the results of the numerically exact 
hierarchical equations of motion (HEOM) approach, we show how the inclusion of vibrational modes 
influences the transfer. 



I. INTRODUCTION 

Since the pioneering work of Franck and Teller [l| , ex- 
citon theory is widely used to describe the optical and 
energy transfer properties of light harvesting systems . 
In recent years particular focus has been drawn to the so- 
called Fenna-Matthews-Olson (FMO) complex (see e.g. 
Refs. 0423) J that promotes energy transfer from the 
main light harvesting complex towards the reaction cen- 
ter in green sulfur bacteria 0. 

The FMO complex consists of three identical subunits, 
called "monomers" , each containing eight [13] bacte- 
riochlorophyll (BChl) molecules (sec Fig. [1]). It is as- 
sumed that BChls 1, 6, and 8 are located at the base- 
plate and receive electronic excitation captured by the 
chlorosomes [13] ■ BChls 3 and 4 are located in the vicin- 
ity of the reaction center. To understand the excitation 
transfer mechanisms, many experimental and theoretical 
investigations, that took advantage of the availability of 
high resolution crystal structure, have been performed. 



FIG. 1: The numbering of the BChls within the FMO 
monomer. In the present work BChl 8 is not taken into ac- 
count. The figure is created using PyMOL [2^ . 
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Supported by theoretical modeling a good overall under- 
standing of the local energies of the chlorophyll molecules 
and their dipole-dipole couplings has been obtained, al- 
though there is still a large variance between values found 
in the literature H^. 

The present interest in the FMO complex has been 
greatly stimulated by the evidence for wavelike energy 
transfer through quantum coherence in photosynthetic 
systems [lo| . Since the electronic coherences arc strongly 
influenced by the coupling of electronic excitation to vi- 
brational modes of the BChl molecules and by the cou- 
pling to the protein environment, a detailed modeling of 
this interaction is important. Information on this cou- 
pling has been obtained e.g. from fluorescence line nar- 
rowing d, [2^ or from theoretical simulations [l^ [U, , 
which indicate that the spectral density, describing the 
coupling of the electronic excitation to vibrational modes, 
is quite structured. Such a structured spectral density 
will lead to non-Markovian effects. The influence of these 
non-Markovian effects have been investigated by various 
theoretical methods [30l - l36j . In general these methods 
are difficult to apply to large systems and/or arbitrary 
spectral densities. 

In the present paper we present an efficient approach, 
which allows us to treat structured spectral densities 
with minimal requirements on the computer capabilities. 
All the calculations presented below can be performed 
on a standard PC within a few minutes. The method 
is based on the non-Markovian quantum state diffusion 
(NMQSD) approach [13-111 and allows us to treat the 
whole range from coherent wavelike motion to incoher- 
ent hopping. We have applied this method, within the 
zeroth order functional expansion (ZOFE) approxima- 
tion, to describe the optical and transfer properties of 
molecular aggregates j43l - |45l |. Within the ZOFE approx- 
imation it is possible to derive a non-Markovian master 
equation [40j . which we use in this paper (in the follow- 
ing we will call it ZOFE master equation). This approach 
enables us to calculate the excitation dynamics of fairly 
large systems, e.g. aggregates consisting of more than 
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'^SO chromophorcs. taking structured spectral densities 
for the exciton-phonon coupling into account. 

Although our method contains an approximation, from 
a previous study on absorption spectra of smaller com- 
plexes we expect that for parameters relevant for the 
FMO complex the approximation works well [i^. This 
will be confirmed in the following, where we will compare 
excitation energy transfer, calculated using the ZOFE 
master equation, with results of the so-called hierarchi- 
cal equations of motion (HEOM) approach. The HEOM 
approach has become popular [l^, ^M, I46l - l49| , since it al- 
lows in principle converged results to be obtained. How- 
ever, the numerical effort grows rapidly with increasing 
system size and for spectral densities that deviate from 
the Drude-Lorentz (D-L) form. To overcome these dif- 
ficulties extensive parallclization (e.g. by using graphics 
processing units [l8|) is necessary. 

Since in our method we are not restricted to a sim- 
ple form of the spectral density, we can investigate how 
energy transfer depends on its details. In particular we 
show that the transfer dynamics is not sensitive to cou- 
pling to modes with energies above the largest energy 
difference between the purely electronic eigenstates of 
the FMO complex (which is roughly 500 cm~^ ). Fur- 
thermore, we demonstrate that the transfer dynamics de- 
pends on the details of the local structure of the spectral 
density in the region, where the FMO complex has elec- 
tronic transition energies. 

The paper is organized as follows: In Section |lT] we 
introduce the model used to describe the FMO complex 
and we present the ZOFE master equation that we use 
for the numerical calculations. In Section IIIII we first 
demonstrate that with the ZOFE master equation we are 
able to reproduce calculations obtained with the HEOM 
method. Then we study the infiuence of different parts of 
the spectral density. Finally, in Section |IV] we conclude 
with a summary and an outlook. 



II. MODEL 

In the following we will briefiy introduce the model 
Hamiltonian that we use to describe the FMO complex. 
Further details can be found e.g. in Ref. j45l[50|. 

In the equations below, we set h = 1. For each chro- 
mophore we take two electronic states into account. Since 
we are interested in excitation transfer of a single elec- 
tronic excitation, we restrict ourselves to the one-exciton 
subspace which is spanned by the states 
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where chromophore n is excited electronically and the 
other BChls are in their electronic ground state. In this 
basis the electronic "system" part of the FMO complex 
is 
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where e„ are the electronic energies of the chromophorcs 
and Vnm are the couplings between them. The "environ- 
ment" of vibrational modes is taken in the form 
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and the coupling of electronic excitation to these vibra- 
tions is expressed through 
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Here the k„a, which describe the coupling of the elec- 
tronic excitation on BChl n to the local vibrational mode 
A at this chromophore, are related to the spectral density 
J{u;) by [13 
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which will be taken as a continuous function in the fol- 
lowing. For a thermal initial state pth of the environment 
the environment correlation function a„(t — t') is related 
to the spectral density by the standard expression 
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coth 



2kT 



cos(ajT) — i sin(a;r) (6) 



Note that we assume that the local environments of dif- 
ferent chromophorcs are uncorrelated, which is in accor- 
dance with recent calculations [13. Finally the total 
Hamiltonian in the one-exciton space is given by 



H" = H, 
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In the following we consider electronic excitation trans- 
fer. Thus we will focus on the reduced density operator 
p{t) in the electronic subspace, which is the trace of the 
total density operator ptot{t) of system and environment 
over the vibrations 



p{t) =^ TronvPtot(i)- 



(8) 



Its diagonal elements in the 1 7r„ ) basis describe the time- 
dependent populations of the chromophorcs. 



A. Numerical method 

We calculate the reduced density operator by applying 
a recently developed method, which is based on the non- 
Markovian quantum state diffusion approach (ssl . l4ll . l43l . 
I45I . Within the ZOFE approximation (described in detail 
in Ref. [s^ H^) one can derive a convolutionless master 
equation for the reduced density operator [4^ 
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Here the operator 
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is determined by the auxihary equation 



(11) 
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with the initial condition Oq (s,s) — — |7r„}(7r„| (we 
use the same notation as in Ref. [43 where the subscript 
indicates that the operators are obtained from the ZOFE 
approximation). For the numerical implementation, we 
approximate the environment correlation function as a 
sum of exponentials, i.e. an{T) = gi^njr -^[fj^ jj^ 

general complex prefactors pnj and complex frequencies 
i7„j. Such a form allows us to derive a closed set of 
coupled differential equations 



Pnj I '^n 



(12) 



for the operators OQ^\t) = J2j O^^"" (t) with initial con- 
dition Oq"'''' (t = 0) = 0. In the numerical implemen- 
tation we expand the equations above in the 1 7r„ ) ba- 
sis. To determine the flnj we proceed similarly as in 
Ref. [51I, [S^l- However, in the present work we do not 
use the Matsubara decomposition for the coth appear- 
ing in the correlation function a„(r) (see Eq. ([5])) but 
use the partial fraction decomposition proposed by Croy 
and Saalmann [53[ . The latter has superior convergence 
properties and uses poles in the complex plane (not only 
on the imaginary axis, as is the case for the Matsubara 
decomposition). Since in the present approach one can 
use hundreds of exponentials, in principle it is no prob- 
lem to approximate the environment correlation function 
air) for arbitrary spectral densities J{uj) and tempera- 
tures with high accuracy. 



III. EXCITATION ENERGY TRANSFER 

A. Comparison with HEOM calculations 

First we will compare the energy transfer calculated 
with the ZOFE master equation approach with the 
HEOM results of Ref. [l^. For the comparison we re- 
strict ourselves to the same model used in Ref. [IH where 
only a single monomer unit of the FMO trimer was con- 
sidered and the recently discovered eighth BChl was not 
taken into account (see Hamiltonian in Table |l|) . 
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-53.5 
320 



6.7 

0.7 

-2.2 
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-13.7 
11.8 
-9.6 
-17.0 
81.1 
630 



-9.9 
4.3 
6.0 
-63.3 
-1.3 
39.7 
440 



TABLE I: 

FMO monomer in cm~^. The upper triangle (couplings) is 
the same as the lower triangle (not shown here). From the 
energies on the diagonal an offset of 12000 cm~^ is subtracted 
(for the calculation of the transfer only the energy differences 
are relevant). The Hamiltonian is that used in Ref. [S^] with 
the energies and couplings from Ref. [l^ . 
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FIG. 2: The different spectral densities considered, (a) Solid 
line: D-L spectral density used in Ref. [H]. Triangles: as 
solid line, but with faster drop-off above 500 cm~^. Dots: 
same as solid line, but with extra peak at 1600 cm~^. (b)- 
(d) Dots: as solid line in (a). Solid line: fictitious structured 
spectral densities, where one peak at a time is enhanced, while 
the effective reorganization energy Aofi is kept constant and 
equal to that of the dotted curve. The vertical lines indicate 
the values of all transition energies between the electronic 
eigenenergies of the FMO monomer unit. 



For the comparison we have taken the transition ener- 
gies and couplings from [l^. In Ref. 32 1 a Drude-Lorentz 
(D-L) spectral density was used to describe the coupling 
to the environment. For the decomposition of the envi- 
ronment correlation function as a sum of exponentials, 
described above, we represent the spectral density by 
a sum of anti-symmetrized Lorentzians. With 10 anti- 
symmetrized Lorentzians we could fit the D-L spectral 
density of Ref. [IH perfectly in the relevant energy range. 
The resulting spectral density is shown as the solid curve 
in Fig.[i;a). 

The calculated transfer resulting from this spectral 
density is shown in Fig. [3] The curves show the 
time-dependent excitation probabilities of the individual 
BChls — the dashed lines are the HEOM results taken 
from Ref. |32| and the solid lines show the results calcu- 
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FIG. 3: Comparison of the excitation transfer taken from 
Ref. [s^] calculated with HEOM approach (dashed lines) and 
that obtained from the ZOFE master equation (solid lines). 
The used spectral density is the solid curve of Fig. EJa). (a) 
Excitation probability for BChl 1 to 4 over time at tempera- 
ture 77K, with initial excitation only on BChl 1. (b) Excita- 
tion probability for BChl 3 to 6, 77K, and initial excitation 
only on BChl 6. (c) As (a), but for 300K. (d) As (b), but for 
300K. 



lated vifith the ZOFE master equation (see Eq. 1^). In 
the upper panels the transfer is calculated for a tempera- 
ture of 77K and in the lower panels for 300K. Initially all 
excitation is either on BChl 1 (left panels) or on BChl 6 
(right panels). In all these four cases there is quite good 
agreement between the HEOM and the ZOFE results, 
showing that the ZOFE approximation preserves all the 
detailed features of the transfer, e.g. the strong oscilla- 
tions at 77K. 



First let us consider the influence of the high energy 
part of the spectral density. As an example we consider 
a spectral density (dotted curve in Fig. EJa)), which co- 
incides with the D-L spectral density of Ref. [s^ for 
low energies, but has a large additional peak centered at 
about 1600 cm~^, resulting in strong additional coupling 
to high energy modes. We have calculated the transfer 
through the FMO monomer for this modified spectral 
density using the ZOFE master equation. The result- 
ing transfer is exactly the same as that calculated by the 
ZOFE method for the original spectral density of Ref. [32| 
shown in Fig. [31 This shows that despite contributing 
significantly to the reorganization energy, the high en- 
ergy part above 500 cm~^ (up to 500 cm~^ the spectral 
densities are the same) is of no relevance for the system 
dynamics. The reason for that lies in the fact that the 
system energies, i.e. all the transition energies between 
the electronic eigenenergies of the FMO monomer (indi- 
cated by the vertical lines in Fig. [2|), are in the energy 
region below 500 cm^^. That means that the detuning 
between these system energies and the high energy part 
of the spectral density is so large (compared to the cou- 
pling strength) that the high energy modes do not couple 
to the system. However, the reorganization energy does 
not take this fact into account and thus gives in this case 
a wrong indication for the strength of the coupling to the 
modes. 

Similarly, a faster drop-off, as given by the spectral 
density shown by the curve (triangles) in Fig. [2ja), does 
not affect the transfer dynamics, since in the region below 
500 cm~^ it also agrees with the original spectral density. 



C. Influence of low-energy modes 



B. Influence of high-energy modes 

Since the ZOFE method enables us to consider an al- 
most arbitrary form of the spectral density, we can now 
investigate the infiucncc of different contributions to the 
spectral density. 

Often, the so-called reorganization energy 

A = / duj J{uj)/uj (13) 

is used as a global measure for the coupling strength to 
the (environmental) vibrational modes depending on the 
spectral density J{uj) Then strong coupling (com- 
pared to the electronic interactions Vnm) is identified as 
A ^ Vnm and weak coupling as A ^ Kim- 

In the following we will investigate the influence of con- 
tributions to the spectral density in different energy re- 
gions and we will show that in many cases the reorganiza- 
tion energy is not a reasonable measure for the coupling 
strength. 



The spectral densities we considered so far have a 
smooth form without any prominent structure. However, 
a more realistic spectral density consists of several dis- 
tinct peaks, that embody the strong coupling to the intra- 
molecular modes of the BChls, and is taken e.g . from 
measured fluorescence line narrowing spectr a |l6l . l29l| or 
from atomistic simulations (see e.g. Ref. [Ill, |28||). 

In the following we will investigate the influence of such 
a structured spectral density and the influence of how the 
coupling strength is shared among the peaks at differ- 
ent energies. To this end we consider a spectral density 
consisting of four distinct peaks. To obtain information 
about the importance of the individual peaks we vary 
their height in a range of parameters where from previ- 
ous investigations [43 we are confident that the ZOFE 
approximation gives reliable results. 

To distinguish the influence of the structure from that 
of the overall coupling strength, we want to roughly keep 
the total coupling strength in the relevant energy region 
(below ~ 550 cm~^) constant. To this end we first intro- 
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FIG. 4: Influence of the structure of the spectral density: 
Transfer calculated with ZOFE master equation, as in Fig. O 
but for the three different spectral densities in Fig. [Hb)-(d) 
and the D-L spectral density in Fig. [2lja). 

duce an effective reorganization energy 

Aoff = I duj J{uj)/uj, (14) 

where we take i^min = cm"^ and iSmax = 550 cni~^ 
as the relevant energy range, that contains all electronic 
transition ener gies of the FMO monomer (see the vertical 
lines in Fig. [2|)|58j. 

For the D-L spectral density in Fig.[2ja) we get a value 
of Aeff = 31 cm~^ for the effective reorganization energy. 

To investigate the influence of the structure of the spec- 
tral density, we consider three different spectral densities, 
that all have the same effective reorganization energy as 
the D-L spectral density. But in contrast to the D-L spec- 
tral density, these spectral densities consist of four promi- 
nent peaks, whose intensities are varied. These spectral 
densities are shown as the solid curves in Fig. [2l^b)-(d). 
From Fig. [2jb) to (d), one of the peaks at a time is en- 
hanced. For comparison again the D-L spectral density 
is shown in each figure as the dotted curve. Note, that 
due to the division by w in the definition of the reorga- 
nization energy (see Eq. (HH)), to keep the value of Aeff 
constant, the peak intensity of the high-energy peaks be- 
comes much larger than that of the low-energy peaks, 
when they are enhanced. 

To investigate the influence of these changes in the 
spectral density, we calculated the excitation transfer 
through the FMO monomer using the ZOFE master 
equation. The resulting time-dependent excitation prob- 
abilities are shown in Fig. 2] in the same way as in 
Fig. [3l but now for the three different spectral densities 
of Fig. [Hb)-(d) and the D-L spectral density. 

One sees that even though the overall coupling 
strength in the relevant energy region is kept constant, 
the transfer changes clearly when the structure of the 
spectral density is varied. In Fig. SJa) for instance, the 



population on BChl 3 after Ips changes almost by up to 
50% depending on the structure of the spectral density. 

Obviously, not only the overall coupling strength in the 
relevant energy region is important, but also the local 
distribution of the coupling strength. 

When the temperature is increased (lower panels), the 
spreading of the curves decreases, showing that the struc- 
ture of the spectral density becomes less important at 
higher temperatures. 



IV. SUMMARY AND OUTLOOK 

In the present paper we have discussed excitation en- 
ergy transfer in the FMO complex, emphasizing the role 
of the coupling of the excitation to vibrational modes of 
the BChl molecules. 

a. ZOFE master equation: To handle the coupling 
of the electronic excitation to vibrational modes the 
numerical calculations have been performed using the 
ZOFE master equation, utilizing the ZOFE approxima- 
tion, which allows to calculate the transport for the FMO 
complex within a few minutes on a standard PC. We 
validated the applicability of the ZOFE approximation 
by comparing with the exact HEOM results of Ref. [s^l , 
where a Drude-Lorentz spectral density has been used. 

b. Reorganization energy: Since we are not re- 
stricted to a certain spectral density (SD), we investi- 
gated the influence of the parts of the SD in different 
energy regions. As expected, the high energy part of the 
SD (above ^ 500 cm~^) does not influence the trans- 
port properties. In this context it is worth noting that 
the so-called reorganization energy, which is often used 
as a measure for the strength of the coupling to the en- 
vironment, is often not a reasonable measure, since it 
takes into account the SD at all energies. In particular 
SDs, that have the same reorganization energy but, in 
the relevant energy region, have different shapes, affect 
the system dynamics differently. 

c. Structured spectral density: Being aware of this 
we used an "effective reorganization energy" (ERE), 
which is a local measure in the energy region of sys- 
tem transitions jsoj , to investigate the dependence of the 
transfer on the local structure of the SD. We find that 
even when the structure is changed only slightly and the 
ERE is kept constant, we observe a change in the transfer 
to BChl 3 up to 50%. This influence of the structure of 
the SD decreases when temperature is increased. 

d. Electronic energies and couplings: Note that the 
transfer of course depends also strongly on the precise 
values of the electronic couplings Vnm between the BChls 
and their transition energies. For example we found that 
the final population of BChl 3 increases by roughly 50% 
when the electronic couplings are increased by 25%. Thus 
care has to be taken not to overestimate theoretical re- 
sults obtained for a certain set of parameters. 

e. Outlook: Since measured spectral densities show 
many peaks in the energy region where electronic tran- 
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sition energies are located [291 ] . in further work the in- 
terplay between the electronic Hamiltonian and highly 
structured SDs should be studied in more detail. 

For such systematic studies the NMQSD (and accord- 
ingly the ZOFE master equation derived from it) seems 
to be well suited, due to the short calculation times. An- 
other big advantage of the computational efficiency of the 
approach is that larger systems arc tractable now - the 
full trimer with 24 chromophores and a structured SD 
can be calculated on a standard PC within a few hours. 
One goal of the present study was to further investigate 
the rather abstract ZOFE approximation [H, H^. The 
quite good agreement with the exact HEOM results of 
Ref. [321 makes us confident that it is worth to further 
explore the NMQSD method. 

The NMQSD approach is also suitable to describe 
other light harvesting complexes and molecular aggre- 
gates. While in the case of the FMO complex the cou- 
pling to high energy vibrational modes with energies 
larger than 1000 cm^^ does not influence the energy 
transfer markedly, in molecular aggregates of organic 
dyes these modes have a strong impact on the electronic 
excitation dynamics [s^ . Furthermore they have ener- 
gies comparable to the dipole-dipole interaction which 
influences strongly e.g. the aggregate absorption spectra 
[55| . In previous investigations we have found that typ- 



ical features of molecular aggregate spectra, such as the 
exchange narrowing of the J-band and the broad H-band, 
are reproduced by the NMQSD calculations 

The NMQSD method is also readily applicable to treat 
optical spectroscopy. In particular it should be well 
suited to investigate the influence of the coupling to vi- 
brational modes on the signals observed in 2D experi- 
ments. Such a study might be revealing, since the "power 
spectrum" obtained from 2D measurements in Ref. 
shows strong resemblance to the fluorescence line nar- 
rowing spectrum of Wcndling et. al [29j . 

To obtain further information how individual vibra- 
tional modes influence the transport it might be useful 
to simulate the transport using networks consisting of 
quantum [s^ or classical oscillators [STj . 
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